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Estimation of Gauss-Markov Random Fields 
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Abstract 

We present telescoping recursive representations for both continuous and discrete indexed noncausal 
Gauss-Markov random fields. Our recursions start at the boundary (for example, a hypersurface in R d , 
d > 1) and telescope inwards. Under appropriate conditions, the recursions for the random field are 
differential/difference representations driven by white noise, for which we can use standard recursive 
estimation algorithms, such as the Kalman-Bucy filter and the Rauch-Tung-Striebel smoother. 

Index Terms 

Random Fields, Gauss-Markov Random Fields, Reciprocal Processes, Gauss-Markov Random Pro- 
cesses, Kalman Filter, Rauch-Tung-Striebel Smoother, Recursive Estimation, Recursive Processing, Tele- 
scoping Representation 

I. Introduction 

We consider the problem of deriving recursive representations for spatially distributed signals, such 
as temperature in materials, concentration of components in process control, intensity of images, density 
of a gas in a room, stress level of different locations in a structure, or pollutant concentration in a lake 
CD-EL These signals are often modeled using random fields, which are random signals indexed over M. d 
or 7L d , for d > 2. For random processes, which are indexed over M, recursive algorithms are recovered 
by assuming causality. In particular, for Markov random processes, the future states depend only on the 
present state given both the past and present states. When modeling spatial distributions by random fields, 
it is more appropriate to assume noncausality as opposed to causality. This leads to noncausafl Markov 
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'When referring to causal or noncausal Markov random fields or processes, we really mean they admit recursive or nonrecursive 
representations. 
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(a) Causal Structure (b) Noncausal Structure 

Fig. 1. Causal and Noncausal models for random fields 



random fields (MRFs): the field inside a domain is independent of the field outside the domain given the 
field on (or near) the domain boundary. The need for recursive algorithms for noncausal MRFs arises to 
reduce the increased computational complexity due to the noncausality and the multidimensionality of the 
index set. The assumption of noncausality presents problems in developing recursive algorithms such as 
the Kalman-Bucy filter for noncausal MRFs. Instead, to derive recursive algorithms, many authors make 
causal approximations to random fields over M 2 or 1?, see Rl- lfTTTl . An example of a random field with 



causal structure is shown in Fig. |l(a)| It is assumed that the site indicated by 'o' depends on the neighbors 
indicated by 'x'. Such fields do not capture fully the spatial dependence, as for example, when the field 
at a spatial location depends on its neighbors. More appropriate representations are noncausal models, an 
example of which is the nearest neighbor model shown in Fig. |l(b)| In lfT2l . the authors derive recursive 
estimation equations for nearest neighbor models over 1? by stacking two rows (or columns) at a time 
of the lattice into one vector and thus converting the two-dimensional (2-D) estimation problem into a 
one-dimensional (1-D) estimation problem with state of dimension 2n, for annxn lattice. However, the 
algorithm in lfT2l is restricted to nearest neighbor models with boundary conditions being local, i.e., they 
involve only neighboring points along the boundary. In |fT3ll , the authors derive a recursive representation 
for general noncausal Gauss-Markov random fields (GMRFs) over Z 2 by stacking the field in each row 
(or column) and factoring the field covariance to get 1-D state-space models. However, the models in 
|[T3l are only valid when the boundary conditions are assumed to be zero. Further, since we can not stack 
columns or rows over a continuous index space, it is not clear how the methods of |[T2l and |[T3l can be 
extended to derive recursive representations for noncausal GMRFs over M. d for d>2. 

For noncausal isotropic GMRFs over M 2 , the authors in Ifl4ll derived recursive representations, and 
subsequently recursive estimators, by transforming the 2-D problem into a countably infinite number 
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of 1-D problems. This transformation was possible because of the isotropy assumption since isotropic 
fields over M 2 , when expanded in a Fourier series in terms of the polar coordinate angle, the Fourier 
coefficient processes of different orders are uncorrected Ifl4ll . In this way, the authors derived recursive 
representations for the Fourier coefficient process. The recursions in lfl4ll are with respect to the radius 
when the field is represented in polar coordinate form. The algorithm is an approximate recursive 
estimation algorithm since it requires solving a set of countably infinite number of 1-D estimation 
problems [14]. For random fields with discrete indices, nonrecursive approximate estimation algorithms 
can be found in the literature on estimation of graphical models, e.g., |fT31 . 

In this paper, we present a telescoping recursive representation for general noncausal Gauss-Markov 
random fields defined on a closed continuous index set in M d , d > 2, or on a closed discrete index 
set in Z d , d > 2. The telescoping recursions initiate at the boundary of the field and recurse inwards. 
For example, in Fig. [2ta), for a GMRF defined on a unit disc, we derive telescoping representations 
that recurse radially inwards to the center of the field. For the same field, we derive an equivalent 
representation where the telescoping surfaces are not necessarily symmetric about the center of the disc, 
see Fig. [2jb). Further, the telescoping surfaces, under appropriate conditions, can be arbitrary as shown 
in Fig. |2tc). In general, for a field indexed in M. d , d > 2, the corresponding telescoping surfaces will be 
hypersurfaces in R d . We parametrize the field using two parameters: AG [0, 1] and 9 G 6 C IR d_1 . The 
parameter A indicates the position of the telescoping surface and the set parameterizes the boundary of 
the index set. For example, for the unit disc with recursions as in Fig. 12a), the telescoping surfaces are 
circles, and we can use polar coordinates to parameterize the field: radius A and angle 6 G © = [— tt, 7t]. 
The telescoping surfaces are represented using a homotopy from the boundary of the field to a point 
within the index set (which is not on the boundary). 

The key idea in deriving the telescoping representation is to establish a notion of "time" for Markov 
random fields. We show that the parameter A, which corresponds to the telescoping surface, acts as time. 
In our telescoping representation, we define the state to be the field values at the telescoping surfaces. 
The telescoping recursive representation we derive is a first order representation in the parameter A and 
is driven by Brownian motion. For a certain class of homogeneous isotropic GMRFs over M 2 , for which 
the covariance is a function of the Euclidean distance between points, we show that the driving noise is 
2-D white Gaussian noise. For the Whittle field |[T6ll defined over a unit disc, we show that the driving 
noise is zero and the field is uniquely determined using the boundary conditions. 

Using the telescoping recursive representation, we promptly recover recursive algorithms, such as the 
Kalman-Bucy filter ifTTIl and the Rauch-Tung-Striebel (RTS) smoother |fT8ll For the Kalman-Bucy filter, 
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we sweep the observations over the telescoping surfaces starting at the boundary and recursing inwards. 
For the smoother, we sweep the observations starting from the inside and recursing outwards. Although, 
we use the RTS smoother in this paper, other known smoothing algorithms can be used as well, see 

ED-ED- 

We derive the telescoping representation in an abstract setting over index sets in R d , d > 2. We can 
easily specialize this to index sets over Z d , d > 2. We show an example of this for GMRFs defined 
over a lattice. We see that, unlike the continuous index case that admits many equivalent telescoping 
recursions, the telescoping recursion for discrete index GMRFs is unique. 

The organization of the paper is as follows. Section JI] reviews the theory of GMRFs. Section ITTT1 
introduces the telescoping representation for GMRFs indexed on a unit disc. Section HVl generalizes the 
telescoping representations to arbitrary domains. Section [V] derives recursive estimation algorithms using 
the telescoping representation. Section|VI]derives telescoping recursions for GMRFs with discrete indices. 
Section I VII I summarizes the paper. 

II. Gauss-Markov Random Fields 

A. Continuous Indices 

For a random process x(t), t £ R, the notion of Markovianity corresponds to the assumption that the 
past {x(s) : s < t}, and the future {x(s) : s > t} are conditionally independent given the present x(s). 
Higher order Markov processes can be considered when the past is independent of the future given the 
present and information near the present. The extension of this definition to random fields, i.e., a random 
process indexed over M d for d > 2, was introduced in |[22l . Specifically, a random field x(t), t G T C M. d , 
is Markov if for any smooth surface dG separating T into complementary domains, the field inside is 
independent of the field outside conditioned on the field on (and near) dG. To capture this definition in a 
mathematically precise way, we use the notation introduced in ll23l . On the probability space (SI, J r , T 3 ), 
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Fig. 3. (a) An example of complementary sets on a random field defined on T with boundary dT. (b) Corresponding notion 
of complementary sets for a random process. 

lej^| x(t) G K be a zero mean random field for t G T C R d , where d > 2 and let dT C T be the smooth 
boundary of T. For any set A C T, denote x(A) as 

x{A) = {x(t) : t G A} . (1) 

Let G- C T be an open set with smooth boundary dG and let G + be the complement of G_ U dG in T. 
Together, G_ and G + are called complementary sets. Fig. [3ta) shows an example of the sets G_, G + , 
and on a domain Tel 2 . For e > 0, define the set of points from the boundary dG at a distance 
less than e 

dG e = {t G T : d(t, SG) < e} , (2) 

where dG) is the distance of a point i G T to the set of points dG. On G± and 9G, define the sets 

E X (G±) = a(x(G ± )) (3) 
E x {dG) = f| a(a;(SG e )) , (4) 

£>0 

where cr(A) stands for the c-algebra generated by the set A. If x(t) is a Markov random field, the 
conditional expectation of x(s), s £ G_, given S a: (G_) is the conditional expectation of x(s) given 

i.e., ma, m 

^[x(s)|E a> (G_)]=S[x(a)|S a ,(5G0], G_ . (5) 

Equation (f5]) also holds for Markov random processes for complementary sets defined as in Fig. [3tb). In 
the context of Markov processes, the set G_ in Fig. f5^b) is called the "past", G + is called the "future", 
and dG is called the "present". The equivalent notions of past, present, and future for random fields is 

2 For ease in notation, we assume x(t) £ R, however our results remain valid for x(t) £ K n , when n > 2. 
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clear from the definition of GL, G + , and dG in Fig. Sa). 

In this paper, we assume x(t) is zero mean Gaussian, giving us a Gauss-Markov random field (GMRF), 
so the conditional expectation in Q becomes a linear projection. Following [24], the key assumptions 
we make throughout the paper are as follows. 

Al. We assume the index set T C M. d is a connected^] open set with smooth boundary dT. 
A2. The zero mean GMRF x(t) £ L 2 (Q, F, V), which means that x{t) has finite energy. 
A3. The covariance of x(t) is R(t, s), where t, s G T C M d . The function space of R(t, s) is associated 
with the uniformly strongly elliptic inner product 

<u,v>=(D a u 1 a a>p D p v^ T = J D a u(s)a a ^(s)D^v(s)ds, (6) 

a|<m,|/3|<m ^ 

where a a ^ are bounded, continuous, and infinitely differentiable, a = [a\, • • ■ , ay] is a multi-index 
of order \a\ = a\ + • ■ • + and the operator D a is the partial derivative operator 

D a = D" 1 • • • D a d d , (7) 

where Df* = d a * /dtf' for t = [h, . . . , t d ). 

A4. Since the inner product in ^> is uniformly strongly elliptic, it follows as a consequence of A3 that 
R(t, s) is jointly continuous, and thus x(t) can be modified to have continuous sample paths. We 
assume that this modification is done, so the GMRF x(t) has continuous sample paths. 

Under Assumptions A1-A4, we now review results on GMRFs we use in the paper. 

Weak normal derivatives: Let dG be a boundary separating complementary sets G_ and G + . Whenever 

we refer to normal derivatives, they are to be interpreted in the following weak sense: For every smooth 

fit), 

y(s) = ^-x{s)^ [ f(s)y(s)dl = lim ^- / f(s)x(s + hs)dl , (8) 

where dl is the surface measure on dG and s is the unit vector normal to dG at the point s. 
GMRFs with order m: Throughout the paper, unless mentioned otherwise, we assume that the GMRF 
has order m, which can have multiple different equivalent interpretation: (i) the GMRF x(t) has m — 1 
normal derivatives, defined in the weak sense, for each point t S dG for all possible surfaces dG, 
(ii) the cr-algebra T, x (dG) in (0]), called the germ cr-algebra, contains information about m — 1 normal 
derivatives of the field on the boundary dG |[24l . or (iii) there exists a symmetric and positive strongly 

3 A set is connected if it can not be divided into disjoint nonempty closed set. 
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elliptic differential operator Ct with order 2m such that ll25l 



C t R{t,s) = 5(t - s) 



(9) 



where the differential operator has the form, 



C t u{t)= (-l) lal D a [a Q ^(t)D^u(t))]- 



(10) 



H,|/3|<m 



Prediction: The following theorem, proved in f24l . gives us a closed form expression for the conditional 
expectation in ([5]). 

Theorem 1 (K2M): Let x(t), t G T C R d , be a zero mean GMRF of order m and covariance R(t, s). 
Consider complementary sets G*_ and G + with common boundary dG. For s g' G_, the conditional 
expectation of x(s) given Y< X (G-) is 



where d^ jdn? is the normal derivative, defined in ([8]), <iZ is a surface measure on the boundary dG, and 
the functions bj(s,r), s ^ G_ and r £ dG, are smooth. 

Proof: A detailed proof of Theorem Q] can be found in (24), where the result is proved for the case 
when s G G + . To include the case when s G 9G, we use the fact that R(t,s) is jointly continuous 



Theorem Q] says that for each point outside GL , the conditional expectation given all the points in GL 
depends only the field defined on or near the boundary. This is not surprising since, as stated before, 
E[x(s)\Tj x (G-)] = E[x(s)\T lx (dG)], and we mentioned before that T lx {dG) has information about the 
m — 1 normal derivatives of x(t) on the surface dG. Appendix IA1 shows how the smooth functions bj (s, r) 
can be computed and outlines an example of the computations in the context of a Gauss-Markov process. 
In general, Theorem Q] extends the notion of a Gauss-Markov process of order m (or an autoregressive 
process of order m) to random fields. 

A simple consequence of Theorem Q] is that we get the following characterization for the covariance 
of a GMRF of order m. 

Theorem 2: If t G G_ and s ^ G_, the covariance R(t, s) can be written as, 




m—1 



(11) 



(consequence of A3) and the uniform integrability of the Gaussian measure (see [26]). 




m—1 



(12) 
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where the normal derivative in ([T2l is with respect to the variable r. 

Proof: Since x(s) - E[x(s)\Y> x (G-)] _L x(t) for t G G-, using (fTTI) . we can easily establish ([f2l . 

■ 

Theorem [2] says that the covariance R(s, t) of a GMRF can be written in terms of the covariance of 
the field defined on a boundary dividing s and t. Both Theorems Q] and |2] will be used in deriving the 
telescoping recursive representation. 

B. Discrete Indices 

Discrete index Markov random fields, also known as undirected graphical models, are characterized 
by interactions of an index point with its neighbors. In this paper, we only consider GMRFs defined on 
a lattice To = [0, JV + 1] x [0, M + 1]. An index G To will be called a node. If two nodes are 

neighbors of each other, we represent this relationship by connecting them with an edge. A path is the 
set of distinct nodes visited when hopping from node to a node {12,32) where the hops are only 

along edges. A subset of sites C separates two sites (i\,ji) £ C and (12,32) ^ C if every path from 
to (12,32) contains at least one node in C. Two disjoint sets A,B C T\C are separated by C if 
every pair of sites, one in A and the other in B, are separated by C. 

We denote the discrete index random field by x(i,j) € R. Let M denote the neighborhood structure 
for the random field, then x(i,j) is a GMRF if x(i,j) is independent of x (Tq\{JV U (i,j)}) given x(N) 
for (i,j) G To\dTo, where OTq denotes the boundary nodes of To. An equivalent way to define GMRFs 
is using the global Markov property: 

Theorem 3 (Global Markov property fiZ7$): For a GMRF x(i, 3) for (i, j) G T = [0, N + l] x [0, M + 
1], for all disjoint sets A, B, and C in To, where A and B are non-empty and C separates A and B, 
x(A) is independent of x(B) given x(C). 

For ease in notation and simplicity, we only consider second order neighborhoods denoted by the set 
A/2 such that for node (0, 0): 

M 2 = {(-1, 0), (1, 0), (0, -1), (0, 1), (1, ±1), (±1, 1), (1, 1), (-1, -1)} . (13) 

Examples of higher order neighborhood structures are shown in Fig. [4] A nonrecursive representation, 
derived in |28], for x(i,j) is given as follows: 

ai,jx(i,j) = ^2 Pij 1 ^ - k J + v^J) > (hj) G T o\dT , (14) 
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• 5 4 3 4 5 



Fig. 4. Neighborhood structure from order 1 to 5. 



where Vij is locally correlated noise such that 



E[v(i,j)x(k,l)]=6(i-k)6(j-l) 

( o (k - i,l - j) i N 2 

E[v(iJ)v(i,j)] = { a id k = i,j = l 

Q% '•' 1 (k-i,l-j)etf 2 



Since E[v (i, j)v(k, I)] = E[v(k,l)v(i, j)], we have 



nk—il—j a i—k,k—l 



(15) 



The boundary conditions in (fl~4l) are assumed to be Dirichlet such that x(dTo) is Gaussian with zero 
mean and known covariance. 



III. Telescoping Representation: GMRFs on a Unit Disc 

In this Section, we present the telescoping recursive representation for GMRFs indexed over a domain 
TcM 2 , which is assumed to be a unit disc centered at the origin. The generalization to arbitrary domains 
is presented in Section [iVl To parametrize the GMRF, say x(t) for t G T, we use polar coordinates such 
that x\{9) is defined to be the point 

x x (0) = x((l - A) cos 0, (1 - A) sin 9) , (16) 

where (A,6>) G [0,1] x [— it, tt]. Thus, {xo(9) : 9 G [— vr,7r]} corresponds to the field defined on the 
boundary of the unit disc, denoted as dT. Let dT x denote the set of points in T at a distance 1 — A 
from the center of the field. We call dT x a telescoping surface since the telescoping representations we 
derive recurse these surfaces. The notations introduced so far are shown in Fig. [5] 
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x{dT Xi ) 




Fig. 5. A random field defined on a unit disc. The boundary of the field, i.e., the field values defined on the circle with radius 
1 is denoted by x(dT). The field values at a distance of 1 — A from the center of the field are given by x(dT^). Each point is 
characterized in polar coordinates as x\(ff), where 1 — A is the distance to the center and 9 denotes the angle. 



A. Main Theorem 

Before deriving our main theorem regarding the telescoping representation, we first define some 
notation. Let x(t) G M be a zero mean GMRF defined on a unit disc T C ffi 2 parametrized as x\(8), 
defined in (fl6l) . Let = [— it, tt] and denote the covariance between x Xl (f?i ) and x\ 2 (9 2 ) by R\ u \ 2 (#1, 9 2 ) 
such that 

Rx u \MM = E[x x ,{9i>x 2 ip 2 )] . (17) 

Define C x {9 1 ,9 2 ) and B x [9) as 

C X (9 X ,9 2 ) = lim ILr x (Q u e 2 )- lim ^-R^xiPufa) (18) 




K C A (c?,c?)=0 

where K is any non-zero constant. We will see in (11261 ) that C\{9, 9) is the variance of a random variable 
and hence it is non-negative. Define Fg as the integral transform 

m_1 f d & 

Fe[x\{B)] = Y i lim -E-bjiM, (A, a))— x x (a)da, (20) 
^ J ti^x+ dp, dni 

where bj((/j,, 9), (A, a)) is defined in (fTTT ) and the index (//, f?) in polar coordinates corresponds to the 
point ((1 — jj) cos 9, (1 — /x) sin 9) in Cartesian coordinates. We see that Fg[x x (9)] operates on the surface 
dT x such that it is a linear combination of all normal derivatives of x x (9) up to order m— 1. The normal 
derivative in (l20b is interpreted in the weak sense as defined in ©. We now state the main theorem of 
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the paper. 

Theorem 4 (Telescoping Recursive Representation): For the GMRF parametrized as x\(9), defined 
in ( fl6l ), we have the following stochastic differential equation 

Telescoping Representation: dx\(6) = Fo[x\(9)]dX + B\(0)dw\(6) , (21) 

where dx\(9) = x\ + d\{9) — x\(9) for dX small, Fq is defined in (1201 ). B\(9) is defined in ( fT9l ), and 
w\(9) has the following properties: 

i) The driving noise w\(9) is zero mean Gaussian, almost surely continuous in A, and independent of 
x(dT) (the field on the boundary). 

ii) For all 6 G 0, w (9) = 0. 

iii) For < Ai < Ai < A 2 < A' 2 and 9 1 ,9 2 G G, wy^i) - w Xl (9i) and wy 2 (9 2 ) - wx 2 {9 2 ) are 
independent random variables. 

iv) For 9 G 0, we have 

v) Assuming the set {u G [0, 1] : C u (9,9) = 0} has measure zero for each 9 G 0, for Ai > A 2 and 
#i,6>2 G 0, the random variable w\ 1 (9) — w\ 2 (9) is Gaussian with mean zero and covariance 



e' 



(wxM - wx 2 (9 2 )) 2 



Proof: See Appendix A. ■ 
Theorem |4] says that x\ + d\{9), where dX is small, can be computed using the random field defined on 
the telescoping surface dT x and some random noise. The dependence on the telescoping surface follows 
from Theorem [TJ The main contribution in Theorem [4] is to explicitly compute properties of the driving 
noise w\(9). We now discuss the telescoping representation and highlight its various properties. 

1) Driving noise w\(9): The properties of the driving noise w\(9) in (1211 ) lead to the following 
theorem. 

Theorem 5 (Driving noise w\(0)): For the collection of random variables 

{w x (9) : (X,9) G [0,1] x 6} 

defined in (1211 ). for each fixed 9 G 0, w\(9) is a standard Brownian motion when the set {u G [0, 1] : 
C u (9,9) = 0} has measure zero for each 9 G 0. 



DRAFT 



12 



Proof: For fixed 9 £ G, to show w\(8) is Brownian motion, we need to establish the following: (i) 
w\(6) is continuous in A, (ii) wq(9) = for all 9 G G, (iii) w\(8) has independent increments, i.e., for 
< X\ < X[ < A2 < A' 2 , w\' i (9) — w\ 1 {9) and w\> 2 (9) — w\ 2 (9) are independent random variables, and 
(iv) for Ai > A2, w\ 1 (6) — w\ 2 (9) ~ AA(0, Ai — A2). The first three points follow from Theorem [4] To 
show the last point, let 9\ = 92 in ( f23T ) and use the computations done in dl34| )-( [T37T ). ■ 
Theorem [5] says that for each fixed 9, w\(9) in (|2TI ) is Brownian motion. This is extremely useful since 
we can use standard Ito calculus to interpret (|2"T1 ). 

2) White noise: A useful interpretation of w\(9) is in terms of white noise. Define a random field 
v\{9) such that 

w x (9)= f u 7 (0)d 7 . (24) 

Using Theorem @1 we can easily establish that v\(0) is a generalized process such that for an appropriate 
function 



C\(9l, (72 

L'<M<7i/'«A^2A|U7 = 

which is equivalent to the expression 



jT *(7)^K(^lK(02)]d 7 = *(A)^p^y, (25) 



%(^(%)1 = ^(Ai - A 2 ) „ C ^ ( ^^ • (26) 
Using the white noise representation, an alternative form of the telescoping representation is given by 

= F e [x x (9)} + B x (9)v x {9) . (27) 

aX 

3) Boundary Conditions: From the form of the integral transform Fq in ( f20t . it is clear that boundary 
conditions for the telescoping representation will be given in terms of the field defined at the boundary 
and its normal derivatives. A general form for the boundary conditions can be given as 

m— 1 



(28) 



y / c k AO, a)— T x (a)da = p k (9) , 9 e @ ,k = 1, . . . ,m , 

where for each k, is a Gaussian process in 9 with mean zero and known covariance. 

4) Integral Form: The representation in (1211) is a symbolic representation for the equation 

xxM = xxM + F e [ Xll {0)W + B^dw^O) , Ai > A 2 . (29) 

Since from Theorem [5J w\(9) is Brownian motion for fixed 9, the last integral in d29l ) is an Ito integral. 
Thus, to recursively synthesize the field, we start with boundary values, given by (l28l) . and generate the 
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field values recursively on the telescoping surfaces dT x for A G (0, 1]. 

5) Comparison to [14]: The telescoping recursive representation differs significantly from the recursive 
representation derived in lfT4"1 . Firstly, the representation in [14] is only valid for isotropic GMRFs and 
does not hold for nonisotropic GMRFs. The telescoping representation we derive holds for arbitrary 
GMRFs. Secondly, the recursive representation in [ 14] was derived on the Fourier series coefficients, 
whereas we derive a representation directly on the field values. 



B. Homogeneous and Isotropic GMRFs 

In this Section, we study homogeneous isotropic random fields over M? whose covariance only depends 
on the Euclidean distance between two points. In general, suppose -R^a^i, #2) is the covariance of a 
homogeneous isotropic random field over a unit disc such that the point (11,61) in polar coordinates 
corresponds to the point ((1 — //)cos0, (1 — //)sin0) in Cartesian coordinates. The Euclidean distance 
between two points (/U,0i) and (A, 62) is given by 

D^ x (0i,e 2 ) = [(1 - /i) 2 + (1 - A) 2 - 2(1 - M )(l - A) cos(0! - 2 )] 1/2 . (30) 

If Rn\(9x, 62) is the covariance of a homogeneous and isotropic GMRF, we have 

R^,x(0i,9 2 ) = r{D fliX (0 1 ,e 2 )) , (31) 

where Y(-) : R — > R is assumed to be differentiable at all points in R. The next Lemma computes 
C\(6i,02) for isotropic and homogeneous GMRFs. 

Lemma 1: For an isotropic and homogeneous GMRF with covariance given by (PTT ). Ca(0i, 02), defined 
in ( fT8l ), is given by 

0i ^ 2 

-2T'(0) 0i = 2 



C x (9 



1,02 



(32) 



Proof: For 9\ 7^ 02, we have 



^-^^(01,02)--! (V tlX (9i,92)) , , , (33) 

where T'(-) is the derivative of the function T(-). Using CU), C A (0i,0 2 ) = when X ^ 2 . 

For 0i = 02, D^ \(9\,92) = \p- — A|, so we have 



^-R,A9i,9 2 ) = -T'(|A - . (34) 

ou I A — p 
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Using CED, Cx(9i,9 2 ) = -2T'(0) when 9 X = 9 2 . ■ 
Using Lemma [T] we have the following theorem regarding the driving noise w x (9) of the telescoping 

representation of an isotropic and homogeneous GMRF. 

Theorem 6 (Homogeneous isotropic GMRFs): For homogeneous isotropic GMRFs, with covariance 

given by (|3TT >. such that T(-) is differentiable at all points in R and Y'(0) < 0, the telescoping 

representation is 



dw x (9) = F e x x (9) + y/-?'(0)dw x (e) . (35) 

For each fixed 9, w\(9) is Brownian motion in A and 

E[w Xl (9 1 )w X2 (9 2 )}=0, Ai ^A 2 ,0i^0 2 (36) 

E[w x (9 1 )w x (9 2 )}=0, 9 X ^9 2 . (37) 



Proof: Since we assume T'(0) < 0, thus B x (9) = C x (9,9) = y/-T'(0), which gives us ([35]). To 
show d36l ) and 071 ), we simply substitute the value of C x (9\,9 2 ), given by (|32l ), in (l22l and use the 
independent increments property of w x (9) given in Theorem 0] ■ 
Example: We now consider an example of a homogeneous and isotropic GMRF where T'(0) = and 
thus the field is uniquely determined by the boundary conditions. Let T(t), t G [0, oo), be such that 

f°° b 

T ( t )= / n ^ h 2^ Mbt)db, (38) 
Jo i 1 + o ) 

where J n (-) is the Bessel function of the first kind of order n l|29l . The derivative of T(t) is given by 

/■oo l2 

r{t) = ~J \TTW Jl{ht)dh > m 

where we use the fact that J' Q {-) = -Ji(-) [29]. Since Ji(0) = 0, T'(0) = and thus B x (9) = in the 
telescoping representation. This means there is no driving noise in the telescoping representation. The 
rest of the parameters of the telescoping representation can be computed using the fact |30l 



(A-l) 2 R(t,s) =S(t-s), (40) 

where R(t,s) corresponds to the covariance associated with T(-) written in Cartesian coordinates and 
A is the Laplacian operator. Since the operator associated with R(t, s) in (l40b has order four, it is clear 
that the GMRF has order two. The field with covariance satisfying (l40l is also commonly referred to as 
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the Whittle field [16]. The telescoping recursive representation will be of the form 



dx x {9) = 




\ d d d ] 

lim —b ((p,9),(X,a))xx{a) + —b 1 ((p,9),(X,a))—x x {a) dadX, 0<A<1, 
t->A+ \_ou ou on 

(41) 



with appropriate boundary conditions defined on the unit circle. 



IV. Telescoping Representation: GMRFs on arbitrary domains 



In the last Section, we presented telescoping recursive representations for random fields defined on a 
unit disc. In this Section, we generalize the telescoping representations to arbitrary domains. Section |TV-A| 
shows how to define telescoping surfaces using the concept of homotopy. Section IIV-BI shows how to 
parametrize arbitrary domains using the homotopy. Section IIV-CI presents the telescoping representation 
for GMRFs defined on arbitrary domains. 

A. Telescoping Surfaces Using Homotopy 

Informally, a homotopy is defined as a continuous deformation from one space to another. Formally, 
given two continuous functions / and g such that f,g:X^y,a homotopy is a continuous function 
h : X X [0, 1] — )■ y such that if x G X, h(x, 0) = f(x) and h(x, 1) = g(x) ED- An example of the use 
of homotopy in neural networks is shown in |[32ll . 

In deriving our telescoping representation for GMRFs on a unit disc in Section [Till we saw that the 
recursions started at the boundary, which was the unit circle, and telescoped inwards on concentric circles 
and ultimately converged to the center of the unit disc. To parametrize these recursions, we can define 
a homotopy from the unit circle to the center of the unit disc. In general, for a domain T C l d with 
smooth boundary dT, the telescoping surfaces can be defined using a homotopy, h : dT x [0, 1] — > c, 
from the boundary dT to a point c G T such that 



PI. {h(t, 0) : t G dT} = dT and {h(t, 1) : t G dT} = c. 

P2. For < A < 1, {h(t, X) : t G dT} C T is the boundary of the region {h(t, p) : (t, p) G dT x (A, 1)}. 
P3. For Ai < A 2 , {h(t,X{),t G dT} c {h{t,p),t G <9T,0 < p < A 2 }. 
P4. \J x {h(t,X),tedT} = T. 

Property 1 says that, for A = 0, we get the boundary dT and for A = 1, we get the point c G T, which 
we choose arbitrarily. Property 2 says that for each A, we want the telescoping surfaces to be in T and 
it should be a boundary of another region. Property 3 restricts the surfaces to be contained within each 
other, and Property 4 says that the homotopy must sweep the whole index set T. 
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Fig. 6. Telescoping Surfaces defined using different homotopies. 



Using the homotopy, for each A, we can define a telescoping surface dT such that 

dT x = {h(9,X) : 9 G dT}, (42) 

where dT is the boundary of the field. As an example, we consider defining different telescoping surfaces 
for the field defined on a unit disc. The boundary of the unit disc can be parametrized by the set of 
points 

dT = {(cos 9, sin 9), 9 G [— 7r, 7r]} . (43) 

We consider four different kinds of telescoping surfaces: 

a) The telescoping surfaces in Fig (6^a) are generated using the homotopy 

h{{cos 9, sin 9), A) = ((1 - A) cos 9, (1 - A) sin 9) . (44) 

b) The telescoping surfaces in Fig Ob) can be generated by the homotopy 

/i((cos0,sin0),A) = ((1 - A)(cos 9 - c x ) + ci,(l- A) (sin 6>-c 2 )+c 2 ), (45) 
= ((1 - A)cos(9 + ci - (1 - A)ci,(l - A)cos6» + c 2 - (1 - A)c 2 ) , (46) 

where (ci,c 2 ) is inside the unit disc, i.e., cf + c| < 1. For the homotopy in (l44l) . each telescoping 
surface is centered about the origin, whereas the telescoping surfaces in (l46l ) are centered about the 
point (ci - (1 - A)ci, c 2 - (1 - A)c 2 ). 

c) In Fig [6£a)-(b), the telescoping surfaces are circles, however, we can also have other shapes for the 
telescoping surface. Fig 0c) shows an example in which the telescoping surface is an ellipse, which 
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we generate using the homotopy 

7i((cos 6, sin#), A) = (a\cos6,b\sm8) (47) 

where and b\ are continuous functions chosen in such a way that P1-P4 are satisfied for h. In 
Fig |6tc), we choose a\ = A and b\ = A 2 , 
d) Another example of a set of telescoping surfaces is shown in Fig Hd). From here, we notice that two 
telescoping surfaces may have common points. 

Apart from the telescoping surfaces for a unit disc shown in Fig [Sa)-(d), we can define many more 
telescoping surfaces. The basic idea in obtaining these surfaces, which is compactly captured by defining 
a homotopy, is to continuously deform the boundary of the index set until we converge to a point within 
the index set. In the next Section, we provide a characterization of continuous index sets in M d for which 
we can easily find telescoping surfaces by simply scaling and translating the points on the boundary. 

B. Generating Similar Telescoping Surfaces 

From Section IIV-AI it is clear that, for a given domain, many different telescoping surfaces can be 
obtained by defining different homotopies. In this Section, we identify domains on which we can easily 
generate a set of telescoping surfaces, which we call similar telescoping surfaces. 

Definition 1 (Similar Telescoping Surfaces): Two telescoping surfaces are similar if there exists an 
affine map between them, i.e., we can map one to another by scaling and translating of the coordinates. 
A set of telescoping surfaces are similar if each pair of telescoping surfaces in the set are similar. 

As an example, the set of telescoping surfaces in Fig Oa)-(b) are similar since all the telescoping 
surfaces are circles. On the other hand, the telescoping surfaces in Fig[6]c)-(d) are not similar since each 
telescoping surfaces has a different shape. The following theorem shows that, for certain index sets, we 
can always find a set of similar telescoping surfaces. 

Theorem 7: For a domain T G M, d with boundary dT if there exists a point c G T such that, for all 
t € T U dT and A G [0, 1], (1 — X)t + Ac G T, we can generate similar telescoping surfaces using the 
homotopy 

h(0,\) = (1- A)0 + Ac, 9edT. (48) 

Proof: Given the homotopy in (l48l . the telescoping surfaces are given by dT x = {h(6, A) : 9 G dT}. 
Using d48]), it is clear that dT° = dT and dT 1 = c. Given the assumption, we have that dT x C T for 
< A < 1. Since the distance of each point on dT x to the point c is (1 — A)||# — c||, it is clear that, for 
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(a) (b) (c) (d) 

Fig. 7. Telescoping Surfaces defined using different homotopies. 

Ai < A 2 , dT Xl C {<9T M : < < A2}. This shows that the homotopy in (|48]> defines a valid telescoping 
surface. The set of telescoping surfaces is similar since we are only scaling and translating the boundary 
dT. M 
Examples of similar telescoping surfaces generated using the homotopy in (|48T ) are shown in Fig [T^a) 
and FigHc). Choosing an appropriate c is important to generate similar telescoping surfaces. For example, 
Fig Gib) shows an example where telescoping surfaces are generated using (|48T ). It is clear that these 
surfaces do not satisfy the desired properties of telescoping surfaces. Fig [T^d) shows an example of an 
index set for which similar telescoping surfaces do not exist since there exists no point c for which 
(1 - \)t + Ac G T for all A G [0, 1] and t G T U dT. 

C. Telescoping Representations 

We now generalize the telescoping representation to GMRFs defined on arbitrary domains. Let x(t) 
be a zero mean GMRF, where t G T C M. d such that the smooth boundary of T is dT. Define a set of 
telescoping surfaces dT x constructed by defining a homotopy h(6, A), where 8 G dT and A G [0, 1]. We 
parametrize the GMRF x(t) as x\(0) such that 

x x (9) = x(h(e,X)). (49) 

Denote G = dT and define C x {0i,9 2 ), B x {9), and F e by dJ), C[9]>, and (US), respectively. Although 
the initial definition for these values was for = [— tt, tt] and x\{6) parametrized in polar coordinates, 
assume the definitions in (fT8T ), ( fT9l , and (f20l > are in terms of the parameters defined in this Section. The 
normal derivatives in the definition of Fg for a point x\(9) will be computed in the direction normal to 
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the telescoping surface dT x at the point h(9,X). The telescoping representation is given by 

dx x {9) = F e [x x (e)]dX + B x {9)dw x {9) (50) 

where the w x (9) = w(h(8, A)) is the driving noise with the same properties as outlined in Theorem [4] It 
is clear from d50l ), that the recursions for the GMRF initiate at the boundary and recurse inwards along 
the telescoping surfaces defined using the homotopy h(9, A). Thus, the recursions are effectively captured 
by the parameter A. 

V. Recursive Estimation of GMRFs 

Using the telescoping representation, we now derive recursive equations for estimating GMRFs. Let 
x(t) be the zero mean GMRF defined on an index set T C M. d with smooth boundary dT. Assume the 
parametrization x x {9) = x(h(9, A)), where h(9, A) is an appropriate homotopy and 9 E © = dT. The 
corresponding telescoping representation is given in (l50l . 

Consider the observations, written in parametric form, as 

dyx(9) = G x (9)x x (9)d\ + D x (9)dn x {9) , < A < 1 , (51) 

where G x {9) and D x (9) are known functions with D x {9) ^ 0, yo(9) = for all 9 G 0, n x (9) is standard 
Brownian motion for each fixed 9 such that 

E[n Xl (9 1 )n X2 (9 2 )}=0, \ 1 ^\ 2 ,9 1 ^ 9 2 , (52) 

and n x {9) is independent GMRF x x {9). 

We consider the filtering and smoothing problem for GMRFs. For random fields, because of the 
multidimensional index set, it is not clear how to define the filtered estimate. For Markov processes, 
the filtered estimate sweeps the data in a causal manner, because the process itself admits a causal 
representation. To define the filtered estimate for GMRFs, we sweep the observations over the telescoping 
surfaces defined in the telescoping recursive representation in ( f50b - Define the filtered estimate x x \ x {9), 



error x x \ x {9), and error covariance S\(a, /3) such that 

z A)A (0) ^ E [x x {9)\o{y^9) , < n < A, 9 e 6}] (53) 

xx\x(0) = xx\x(0)-x xlx (9) (54) 

S x (a,P)±E\x x \ x (a)x X \ X (P)]. (55) 
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The set {y^(9) , < \i < A, 9 G ©} consists of the region between the boundary of the field, dT, and the 
surface dT x . A stochastic differential equation for the filtered estimate x x \ x (9) is given in the following 
theorem. 

Theorem 8 (Recursive Filtering of GMRFs): For the GMRF x\(6) with observations y\{6), a stochas- 
tic differential equation for the filtered estimate x x \ x (9), defined in ( f53T >. is given as follows: 

dx X \ X (9) = F e [x x{x (6)]d\ + K e [de x (e)} , (56) 

where e x (6) is the innovation field such that 

D x {9)de x {9) = dy x (9) - G x (9)x xlx (9)d\ , (57) 
Fq is the integral transform defined in (|2Qb and Kg is an integral transform such that 

K e [de x (e)}= [ Q^S x (a,9)de x (a)da, (58) 

where S x (a,9) satisfies the equation 

^S x (a,9) = F a [S x (a,9)}+F e [S x (a,9)}+C x (9,a)- [ 3^5 A (a, P)S X (9, P)d/3 . (59) 
d\ J D{{fi) 

Proof: See Appendix O ■ 
We show in Lemma |2] (Appendix B) that e x (9) is Brownian motion. Thus, d56l) can be interpreted using 
Ito calculus. Since we do not observe the field on the boundary, we assume that the boundary conditions 
in (l56l ) are zero such that: 

-^jx x]x (e) = o, j = i ...,m-i, ee@. m 

The boundary equations for the partial differential equation associated with the filtered error covariance 
is computed using the covariance of the field at the boundary such that 

d j d j 

—S {a,9) = -^Ro, (a,0), a,9e@. (61) 



The filtering equation in (1561 ) is similar to the Kalman-Bucy filtering equations derived for Gauss- 
Markov processes. The differences arise because of the telescoping surfaces. Using (l56l ). let be a 
single point instead of [— tt,tt}. In this case, the integrals in (|2Qb and (|58l ) disappear and we easily 
recover the Kalman-Bucy filter for Gauss-Markov processes. 

Using the filtered estimates, we now derive equations for smoothing GMRFs. Define the smoothed 
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estimate x x \ T (9), error x x \t{9), and error covariance SW(a, /3) as follows: 

x xlT (6) ^ E[x xlT \a{y(T)}] (62) 

x A |tW = ^aW-Ja|tW (63) 
S MT (a,P)=E[x MT (a)x MT {P)]. (64) 

A recursive smoother for GMRFs, similar to the Rauch-Tung-Striebel (RTS) smoother, is given as follows. 

Theorem 9 (Recursive Smoothing for GMRFs): For the GMRF x x (9), assuming S\(9,9) > 0, the 
smoothed estimate is the solution to the following stochastic differential equation: 

dx x]T (9) = F e [x x]T (9)]dX + g^[£ A[T (0) - x x]x (9)} , 1 > A > , (65) 

where x x \ x (9) is calculated using Theorem[8]and the smoother error covariance is a solution to the partial 
differential equation, 

9S MT (a,0) C x (a,a)S x (a,d) C x (9,9)S x (a,9) 

oT = Fe^AiTia, (9) + F a S x \ T (a, 9) + C A (a, 9) — r , (66) 

dA oa (a, «) &\{9,9) 

where 

Fp = F(, + C x (P,P)/Sx(J3,P). (67) 

Proof: See Appendix |Pl ■ 
The equations in Theorem [9] are similar to the Rauch-Tung-Striebel smoothing equations for Gauss- 
Markov processes [18]. Other smoothing equations for Gauss-Markov processes can be extended to apply 
to GMRFs. 

VI. Telescoping Representations of GMRFS: Discrete Indices 

We now describe the telescoping representation for GMRFs when the index set is discrete. For simplic- 
ity, we restrict the presentation to GMRFs with order two. Let {x(i, j) G R : (i, j) G T\ = [1, N] x [1, M]} 
be the GMRF. Stack each row of the field and form an NM x 1 vector x. In the representation (fT4l . 
stack the noise field v(i,j) row wise into an NM x 1 vector v. The boundary values are indexed in a 
clockwise manner starting at the upper leftmost node into a 2(N + 1) + 2(M + 1) x 1 vector x& 0. A 

4 The ordering does not matter as long as the ordering is known. 
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matrix equivalent of (fl4l) is given as 

Ax = AbX-b + v , (68) 

where A is an NM x NM block tridiagonal matrix with block size N x N, A is an iVA/ x 2(iV + 
1) + 2(M + 1) sparse matrix corresponding to the interaction of the nodes in x with the boundary nodes. 
The matrices A and Af, can be evaluated from the nonrecursive equation given in (fT4l ). Further, we have 
the following relationships, 

£[xv T ] = I and E[vv T ] = A . (69) 

Equation (l68l is an extension of the matrix representation given in |[T3l for the case when x& = 0, i.e., 
boundary conditions are zero. For more properties about the structure of the matrix A, we refer to lfi3l . 

Let T k = [k, N + 1 - k] x [k, M + 1 - k] and let dT k be the boundary nodes of the index set T k 
ordered in a clockwise direction. For example, x(dTo) = x&. Define z k such that 

z k 4 x (dT k ) = x(T k \T k+1 ) , (70) 

where k = 0, 1, . . . , [min(M, N) /2] . Define r such that 

t = \mm(M,N)/2] . (71) 

Each z k will be of variable size, and let Mf; be the size of z k , i.e., z k is a vector of dimension x 1. 

As an example, consider the random vectors defined on the 5x5 lattice in Fig. [8] The random vector zq 
consists of the boundary points of the original 5x5 lattice, z\ is the boundary points left after removing 
Zq, and Z2 is the boundary point left after removing both zq and z\. The telescoping nature of zq, z\, 
and Z2 is clear since we start by defining zq on the boundary and telescope inwards to define subsequent 
random vectors. The clockwise ordering of z k is shown by the arrows in Fig. [8] The telescoping recursive 
representation for x is given in the following theorem. 

Theorem 10 (Telescoping Representation for GMRFs): For a GMRF x(i,j), the process z k defined in 
(l70l) is a Gauss-Markov process and thus admits a recursive representation 

Z k = F k Z k _i + w k , k = 1, 2, . . . , T , Zq = x 6 (72) 

where F k is an x Mf_ x matrix and w k is white Gaussian noise independent of z k ^\ such that 

F k = E[z k zZ_ 1 ](E[z k - 1 z%_ 1 ])- 1 (73) 
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.Zl , 















Fig. 8. Telescoping recursions 



Q k = E[w k wl] = E[z k zl] - F k E[z k _ lZ l] . (74) 

Proof: The fact that z k is a Markov process follows from the global Markov property outlined in 
Theorem [3] The recursive representation follows from standard theory of state-space models il9l . ■ 
Both d72l ) and the continuous index telescoping representation are similar since the recursion initiates 
at the boundary and telescopes inwards. For the continuous indices, the recursions were not unique, 
whereas for the discrete index case, the recursions are unique. We outline a fast algorithm for computing 
F k and Q k that does not require knowledge of the covariance of x, just knowledge of the matrices A 
and Ab in d68l 

Define the NM x 1 vector z such that (in Matlab© notation) 

z = [zi;z 2 ; ■■■;zt] ■ (75) 
The random vector z is a permutation of the elements in x such that 

z = Px , (76) 



where P is a permutation matrix, which we know is orthogonal. We can now write ( 1681 ) in terms of z 
by writing x = P T Px = P T z: 

AP T z = A b z + v . (77) 

Multiplying both sides of ( 1771 ) by P, we have 

PAP T z = PA b z + Pv . (78) 

Since E[vv T ] = A, we have E [(Pv)(Pv) T ] = PAP T . This suggests that (|78]> is a matrix based 
representation for the Gauss-Markov process z k . Further, because of the form A and Ab, PAP T and 
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PAb will have the form: 



PAP 7 



M\ -Mf 
-M 2 M° 2 M% 

••• 



PA b = [MJ; 0;... ;0], 



/c-i 



-M 



K 



(79) 



(80) 



where M k is an M k x M|_ 1 matrix, .M° is an x M| matrix, and M k is an M| x Af| +1 matrix. From 
( f69b . PAP T is positive and symmetric, and thus .A/f^ = (.M^ +1 ) T . To find the telescoping representation 
using d78l) . we find the Cholesky factors for PAP T such that 



£i 
-7> 2 C 2 



(81) 



(82) 



• -Vr-l Ct-1 
• • -V T C T 

where the blocks Ck are M k x M k lower triangular matrices, and the blocks Vk are Mf x Mt x matrices. 
Substituting (f8Tb in d78l ) and inverting £ T , we have 



£z = C~ T PAbZQ + C~ T Pv. 



(83) 



Notice that the noise is now white Gaussian since 



E [C~ T Pvv T P T C~ l ] = £- T PAP T C- 1 = I 
If we let V\ = Mj, we can rewrite ( f83l in recursive manner as 

z k = C^VkZk-i + w k , k = 1, 2, . . . , r , 



(84) 



where Ffc = C k l Vk and = C k l C k T . A recursive algorithm for calculating Fk and Qk, which follows 
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from the calculation of the Cholesky factors, is given as follows 1131 : 
Initialization: Q T = {M Q T )~ l , F T = Q T M~ 
For k = t — 1, r — 2, ... ,1 

Qk 1 = Ml - M + k F k+1 
F k = QkM k 

end 

Remark: The telescoping representations we derived shows the causal structure of Gauss-Markov 
random fields indexed over both continuous and discrete domains. Our main result shows the existence 
of a recursive representation for GMRFs on telescoping surfaces that initiate at the boundary of the 
field and recurse inwards towards the center of the field. Just like we derived estimation equations for 
GMRFs with continuous indices, we can use the telescoping representation to derive recursive estimation 
equations for GMRFs with discrete indices. The numerical complexity of estimation will depend on the 
size of the state with maximum size, which for the GMRF is the perimeter of the field captured in the 
state zq. For example, the telescoping representation of a GMRF defined on a yN x y/N lattice with 
non-zero boundary conditions will have a state of maximum size of order 0(\/N). Notice that for both 
continuous and discrete indexed GMRFs, the telescoping representation is not local, i.e., each point in the 
GMRF does not depend on its neighborhood, but depends on the field values defined on a neighboring 
telescoping surface (or F^ is not necessarily sparse). Direct or straightforward implementation of the 
Kalman filter requires 0((\/N) 3 / 2 ) due to a matrix inversion step. However, using fast algorithms and 
appropriate approximations, fast implementation of Kalman filters, see IT331 for an example, can lead to 
0((VN) 2 ), i.e., O(N). 

Now suppose the observations of the GMRF are given by y = Hx. + v , where x G M n is the GMRF, 
if is a diagonal matrix, and v is white Gaussian noise vector such that v ~ A/"(0, R), where R is a 
diagonal matrix. The mmse x is a solution to the linear system, 

[5T 1 +H~ 1 R~ 1 H}x = H T R- l y. (85) 

Since x is a GMRF, it follows from [27] that is sparse, where the non-zero entries in S^ 1 correspond 
to the edges in the grapl^l associated with x. In [35], we use the telescoping representation to derive an 

5 We note that the graphical models considered in this paper are a mixture of undirected and directed graphs, where the 
boundary values connect to nodes in a directed manner. These graphs are examples of chain graphs, see 1341 . and the underlying 
undirected graph can be recovered by moralizing this graph, i.e., converting directed edges into undirected edges and connecting 
edges between all boundary nodes. 
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iterative algorithms for solving (1851 ) using the telescoping representational. Experimental results in |[35l 
suggest that the numerical complexity of the iterative algorithm is O(N), although the exact complexity 
may vary depending on the graphical model. The use of the telescoping representation in deriving the 
iterative algorithm in [35], is to identify computationally tractable local structures using the non-local 
telescoping representation. 

VII. Summary 

We derived a recursive representation for noncausal Gauss-Markov random fields (GMRFs) indexed 
over regions in M. d or Z d , d > 2. We called the recursive representation telescoping since it initiated at the 
boundary of the field and telescoped inwards. Although the equations for the continuous index case were 
derived assuming x(t) is scalar, we can easily generalize the results for x(t) G M n , n > 1. Our recursions 
are on hypersurfaces in M. d , which we call telescoping surfaces. For fields indexed over M. d , we saw that 
the set of telescoping surfaces is not unique and can be represented using a homotopy from the boundary 
of the field to a point within the field (not on the boundary). Using the telescoping representations, we 
were able to recover recursive algorithms for recursive filtering and smoothing. An example of applying 
this to image enhancement of noisy images is shown in [36]. Besides the RTS smoother that we derived, 
other recursive smoothers can be derived using the results in |fT9ll - ll2Tll . We presented results for deriving 
recursive representations for GMRFs on lattices. Extensions of the telescoping representation to arbitrary 
graphical models are presented in |[35l . Using the results in ll35l . we can derive computationally tractable 
estimation algorithms. 

Appendix A 
Computing bj(s,r) in Theorem Q] 

We show how the coefficients bj(s,r) are computed for a GMRF x(t), t G T C M. d . Let G_ and G + 
be complementary sets in T C M. d as shown in Fig. [3] Following [24], define a function h s (t) such that 



h a (t) 



R(s,t) t e G- U dG 

h s {t) : C t h s (t) = , 

(86) 

h s (dG) = R(s,dG) t e G + 
h s (t) e H?(T) 



6 The work in [ 35 1 applies to arbitrary graphical models and the telescoping representations are referred to as block-tree graphs. 
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where C r is defined in (fTOb and H^iT) is the completion of Cq°(T), the set of infinitely differentiable 
function with compact support in T, under the norm Sobolov norm order m. From (O, it is clear that 
h s (r) 7^ R(s,r) when r € G + . Let u(r) G C^°(T) and consider the following steps for computing 

bj(s,r): 



|a],]/3|<m 



^ y D a u{t)a aj3 {t)D p h s {t)dt 

= E / D a u(t)a a> p(t)DPR(s,t)dr+ £ I D a u(t)a a ^(t)D^h s (t)dr (87) 

H,l/3|<m G - |a|,|/3|<rri G + 

= V / bj(s,r)-^u(r)dl + / u(r)L t R(s,t)dt + / u(t)L t h s {t)dt (88) 

= EL b i(^4 u(r)dL (89) 

To get (l87l ). we split the integral integral on the left hand side over G_ and G + . In going from (187T ) 
to ([88]), we use integration by parts and the fact that u{r) £ C^°(T). We get ([89]> using ® and d86l 
Thus, to compute bj(s,r), we first need to find /i s (r) using (l86l ) and then use the steps in (|8"7T)— d89l). 
We now present an example where we compute bj(s,r) for a Gauss-Markov process. 

Example: Let x(t) G R be the Brownian bridge on T = [0, 1] such that 

x(t) = w(t) - tw(l) , (90) 

where w(t) is a standard Brownian motion. Since covariance of w(t) is min(i, s), the covariance of x(t) 
is given by 

s(l-i) t>s 



.R(i,s) = < 



(91) 

t(l-s) t<s 



Using the theory of reciprocal processes, see 11371 . ll38l . it can be shown that the operator £ t is 

L t R(t,s) = - d2R d ^ s) =5(t-s). (92) 

Thus, the inner product associated with R(t, s) is given by 

%1 d , , d 



f x d d 
< u,v >= (Du, Dv) T = \ —u(s)—v(s)ds. 

J ds ds 



(93) 
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Following ([86]), for r < 1 and s G [r, 1], h s {t) = R(s,t) for t G [0,r] and 
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/i s (r) =r(l-s), /i,(l) = 0. 



(94) 
(95) 



We can trivially show that h s (t) is given by 



h a (t) = ^—^-(l-t), t>r. 
1 — r 



We now follow the steps in ([87 



^i^^=H^I^'^ + / Jr^l^ 



u (t)^R( s ,t) 



+ u(t)^h s (t) 



d d 



1 - s 
1 — r 



u(r) . 



Using Theorem [U we can compute i?[x(s)|<r{x(t) : < £ < r}] as 



£[s(s)|a{x(t) : < t < r}] = E[x(s)\x(r)} 



1 - s 
1 — r 



x(r) . 



We note that since x(t) is a Gauss-Markov process, it is known that, [39], 



(96) 



(97) 

(98) 

(99) 
(100) 



(101) 



E[x(s)\x(r)] = R(s, r)R (r,r)x(r) . 



(102) 



Using the expression for R(s,r), we can easily verify that (11011 ) and (11021 ) are equivalent. 

Appendix B 

Proof of Theorem [4j Telescoping Representation 

Let Xx+d\\\(@) denote the conditional expectation of x\ + d\{0) given the u-algebra generated by the 
field {x^(a) : (fi, a) G [0, A] x 0}. From Theorem [TJ we have 

m— 1 „ ,j 

S\+d\\\(0) = / b j ((\ + d\,9),(\,a))—x x (a)da. (103) 
j=0 J ® dnJ 
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It is clear that x x {9) = x\\\{Q). Taking the limit in (1103b as dX — > 0, we have 

m~ 1 



m— 1 ,, 

xaW = X) J_b J ((\,9),(\,a))j-x x (a)da. (104) 



3=0 

Define the error as £, x+dX (9) such that 

£a+c2a($) = X X+dX (9) - X X+dX (9) . (105) 

Adding and subtracting x\{9) in (11051 ) and using (11041 ). we have 

m—l „ .a 

xx+dx{9)-x x {9) = Y J / [b j ((X + dX,9),(X,a))-b j ((X,9),(X,a))}—x x (a)da + ^ x+dX (9). (106) 



e 



3=0 

Assuming dX is small, we can write bj((X + dX, 9), (A, a)) — bj((X, 9), (A, a)) as 

t(\ . w\ m a ^ h rr\ m r\ w ^j((A + rfA, 0), (A, a)) - fej((A,fl), (A, a)) 
6j((A + dA, 0), (A, a)) - bj((X, 9), (A, a)) = — J - dX (107) 



(. lim ^6 i ((M,6>),(A,a))) dA, 



(108) 



where in going from (11071 ) to (1108b . we use the assumption that dA is close to zero. Writing dx\(9) = 
x x+d\(9) — x\(9) and substituting (11081 ) in (11061 ). we get 

dx A (0) = F e [x x (9)]dX + £ A+dA (0) , (109) 

where Fg is given in (l20l . To get the final form of the telescoping representation, we need to characterize 
t\+d\(Q)- To do this, we write £ x+dX (9) as 

Zx +dX (9) = B x (9)dw x {9) = B x (9)[w x+dX (9) - w x (9)} . (110) 

We now prove the properties of w x (9): 
i) Since £\(0) = x x {9) — x\\\{9) and x x \ x {9) = x x {9) by definition, we have 

lim t x+dX (9) = t x (9) = 0,a.s.. 

d\— >0 

Thus, using (1 1 10b . since B x (9) / 0, we have 

lim w x+dX {9) - w x {9) = 0,a.s. (Ill) 

dA-5>0 

lim w x+dX (9) = w x (9),a.s. (112) 

d\— >0 
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Equation (11 121 ) shows that w x (0) is almost surely continuous in A. 

ii) Since the driving noise at the boundary of the field can be captured in the boundary conditions, 
without loss in generality, we can assume that wq(6) = for all 9 G 0. 

iii) For < Ai < X[ < X 2 < A' 2 and 9 1 ,9 2 G 6, let dXi = X[ - X 1 and d\ 2 = X' 2 - A 2 . Consider the 
covariance 

E[i\ 1 +d\ 1 {9l)i\ 2 +d\ 2 (O2] = E [(x Xl+ d Xl (9\) - X Xl+ d Xl (#l)) £,X 2 +dX 2 (02)] (1 13) 

= E[x Xl+dXl (9i)i\ 2+ d\ 2 (02)] - E[x Xl+dXl (9i)(,\ 2 +d\ 2 (02)] (1 14) 
= 0, (115) 

where to go from (II 131 ) to (II 141 ). we use the orthogonality of the error. Using the definition of 
£,\+d\(0) in dHOK we have that (0\) — w Xt (0\) and wy 2 (O2) — w\ 2 (0 2 ) are independent random 
variables. 

iv) We now compute E{£ x +d\(9i)£\+d\(02)Y- 
E [£,x+dx (9i) £,x+dx (02)] 

= E[a x+ dx(0i)(xx + dx(02) - x x+dX[x (9 2 ))] (116) 
= E[(x x+dX (0i) - x A+dA , A (0i))x A+dA (0 2 )] (117) 

m— 1 „ gj 

= Rx+dx,x+dx(0i,02) - E y Q M( A + dA > ^l), ( A > a))^jRx,x + dx(a, 2 )da (118) 
= Rx+dx,x+dx(0i,02) - Rx,x+dx(0i,02) 
+ Rx,x+dx(0i,02)-J2 j_b j ((X + dX,0 l ),(X,a))-^-R XyX+ dx(a,02)da. (119) 



j=0 

Using d 104b , we have 



Rx,x+dx(0i,02) = E[x x (9i)x x+dx (0 2 )] (120) 
= J2 J e b j ((X,9 1 ),(X,a))-^R XMdX (a,e 2 )da. (121) 



j=0 

Substituting (11211 ) in d 1 19b . we have 
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£[&+dA(0l)6+dA(02)] 



R\+d\,\+d\(9i, #2) — R\,\+d\(9i,92 

771 — 1 



^ / [6 i ((A + dA,fli),(A,a))-6 i ((A,0i) > (A,a))]^- J 12 Aj A+dA(a,02)da (122) 
i=o • ye n 



-RA+dA,A+dA(#l, #2) _ R\,\+d\(9l,82) 



(l\ 



dX 



H v dA J iJ^K^ 



^ 777. 1 /1 

lim A (0 l9 2 )_ lim — £ / ft* 



,\a 1 ,a)-r—R\,\{a,6 2 )da | dA 



dA (123) 



(124) 







9 



lim —R x (e u 9 2 )- lim —J? „, A (0i,0 2 ) UA 
jU->-A- ati /*-*A+ a/i / 

= c x (e lf e 2 )d\. 

Thus, for dA small, we have 

£[K+dA(0l) - u/a(0i))(«W 2 ) - tifc(fc))] = f ( f^L dX. 
Since it>o(#) = 0, we can use (1127b to compute E'fu'A (#1)1^ (#2)] as f°H° ws: 



(125) 
(126) 



(127) 



E[(w x (0i) - w {e 1 ))(wx(e 2 ) - w Q (6 2 ))] 



lim E 

iV-*oo 



lim E 

iV-Kx> 



JV+1 iV+l 

K* (*l) - fa)) K* fa) ~ fa)) 

k=0 k=0 
N+l 

$^(™7*fa) ~ ^-ifa))( w 7<cfa) - w^fa)) 

fc=0 



lim > 



C 7fc _i fa, #2) 



(7fc-T*-l) 



n B u (e x )B u {e 2 ) 



du . 



7o = A,7jv+i=0 (128) 

(129) 

(130) 
(131) 



We get d 1291 ) using the orthogonal increments property in (iii). We use d 1271 ) to get d 1 30b - We use 
the definition of the Riemann integrals to go from d 1 301 ) to dl31| ). 
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v) For Ai > A2, the covariance of w\ 1 (9) — w\ 2 {0) is computed as follows: 

E[(w Xl (9) - w X2 (9)) 2 } = E[w 2 Xi (9)] + E[wl (9)} - 2E[w Xl (9)w X2 (9)} (132) 

C u (9i, 9-j) 
B U {9 1 )B U {9 2 ) 



Al + A 2 _ 2 / " 2 n C ^ 2 ) n : du , (133) 



where 







E[w 2 x (9)}= [ X ^Mldu (134) 







% u+ [ ££^ du ,135) 



{ue[0,\]:Cu(9,0)=0} (0) 

du (136) 



{ue[0,\]:C u (8,6)^0} 

= A. (137) 

To go from (TT35T > to (TT361 we use (dU) since C u (9, 9) / 0. To go from (IT36T > to (fT37l ). we use the 
given assumption that the set {u € [0, 1] : C u (9, 9) = 0} has measure zero. 

Appendix C 
Proof of Theorem [8j Recursive Filter 

The steps involved in deriving the recursive filter are the same as deriving the Kalman-Bucy filtering 
equations, see l40l . The only difference is that we need to take into account the dependence of each point 
in the random field on its neighboring telescoping surface (which is captured in the integral transform 
Fq), instead of a neighboring point as we do for Gauss-Markov processes. The steps in deriving the 
recursive filter are summarized as follows. In Step 1, we define the innovation process and show that it 
is Brownian motion and equivalent to the observation space. Using this, we find a relationship between 
the filtered estimate and the innovation, see Lemma [2] In Step 2, we find a representation for the field 
x\{9), see Lemma [5] Using Lemma [2] and Lemma [51 we find a closed form expression for x\\\{9) in 
Step 3. We differentiate this to derive the equation for the filtered estimate in Step 4. Finally, Step 5 
computes the equation for the error covariance. 

Step 1. [Innovations] Define q\(9) such that 

qx(0) = y x (9)- t G^(9)x^(9)d^ . 
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Define the innovation field e\(9) such that 

dex(9) = D^(f) dVx{9) = ^§}^ x{9)dX + dnx{9) ' (138) 

where we have used (I5TT) to get the final expression in (11381) and assume that D\{6) ^ 0. 

Lemma 2: The field e\{9) is Brownian motion for each fixed 9 and E[e\ 1 {9i)e X2 (9 2 )] = when 
Ai ^ A2, 9\ 7^ #2- 

Proo/: Note that E[x x \ x {9)\a{e^{9) : < \i < A}] = since x x \ x (9) _L y M (a) for a G 9, 
< /U < A. Thus, using Corollary 8.4.5 in |40ll , we establish that e\{9) is Brownian motion for each 
fixed 9. Assume Ai > A2 and consider 7 < A2. Then, using the orthogonality of error, _L y 7 (a) 

for 7 < /i, and the fact that n\(9) _L e 7 (a) for 7 < A, we have 

E[{e Xl {0i) -e Aa (e 2 ))e 7 (a)] = jT^ £[5W(0iK(«)] + ^[(«a^i) - n Aa (%))e 7 (a)] 

= 0. (139) 
Now we compute E[de Xl (#i)cte A2 (#2)] for Ai > A2: 



E[de Xl { ai )de X2 {9 2 )] =E 



x Al i Al (0i)dAi + dn Al (0i) de A2 (0 2 ) 



(140) 



EidnxMdexM] (141) 
' £ [dn Al (6>i) {dy X2 (9 2 ) ~ G X2 (9 2 )x X2 \ X2 (9 2 )d\ 2 )] (142) 



1 



■E[dn Xx {0i)dyxM] (143) 



£>A 2 (0 2 ) 

= TTlTr^^ ^ ( G ** &) x ^ ^ d ° 2 + ^ (^)rfn A2 (0 2 ))] (144) 
= E[dn Xl (0 1 )dn X2 (9 2 )] = 0. (145) 

To go from (1140b to (11411 ). we use that x Al | Al (#i) is independent of de X2 {9 2 ), since de X2 {9 2 ) is a linear 
combination on the observations {y(dT s ),s £ [0,A2]}- We get (11421) using the definition of de X2 {9 2 ) in 
(1138b - To go from (11421 ) to (11431 ). we use that dn fll (9i) is independent of x A2 | A2 (#2) for Ai > A2. We get 
(11441 ) using the equation for the observations in (TSTb . To go from (11441 ) to (1 145b . we use the assumption 
that n Al (9\) is independent of the GMRF x x {9). In a similar manner, we can get the result for Ai < A2. 
For Ai / A2 and 9\ / 9 2 , E[e Xl {9\)e X2 (9 2 )] = follows from similar computations as done in 

(DSD-GiB. ■ 

Lemma |2] says that the innovation has the properties as the noise observation n x {9). We now use the 
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innovation to find a closed form expression for the filtered estimate x x \ x (9). 

Lemma 3: The filtered estimate x\\\(8) can be written in terms of the innovation as 

x\\\(0) = / / gx, f ,(0,a)de f ,(a)da (146) 
Jo Je 



d_ 

d\i" 



9\A > «) = — E[x x (e)e^a)] . (147) 



Proof: Using the methods in 11401 or 111911 . we can establish the equivalence between the innovations 
and the observations. Because of this equivalence, we can write the filtered estimate as in (11461) . 
We now compute g XjfM (9,a). We know that 

(xx(6) -x x \\(0)) ± e^(a) , fi<X,a£&. 

Thus, we have 

E[x x (9) efl (a)} = E[x x \ x {9)e^a)\ (148) 

= / / g x . s (0,(3)E[de s ((3)e^a)]d(3 (149) 
Jo Je 

Jo Je Jo 

9\,,(P, (3)5{s - r)5{(3 - a)dsdrdj3 (151) 



o J o Je 

g X} r(6,a)dr. (152) 

o 

To go from ( 11501 ) to dl51| ), we use Lemma|2l Differentiating dl52| ) with respect to /x, we get the expression 

far0 A , M (0,a)inG33. ■ 
Step 2. [Formula for x x (9)] Before deriving a closed form expression for x x {9), we first need the 
following Lemma. 

Lemma 4: For any function ^a^^I; #2) with m — 1 normal derivatives, we have 



f X F gi [^ X}J (0i,e 2 )}d 7 = F dl 
Jo 



'0 

Proof: Using the definition of Fg 1 , we have 



A 



*A, 7 (^1^2)^7 



(153) 



/ F ei [^ Xn (9 1 ,e 2 )}d-/= / Y, / lim — bj((n,9),(\,a)) lim -- ¥ A (a + ha,9 2 )dad 1 
Jo Jo ~q Je m^a+ c/x ft^o cw ' ' 
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m 1 f d (P \ 

V" / li m —bM[i,6),(\,a)) lim -— / 



j=0 



F 6l 



■ A+/iA,7 



(a + /id, #2)^7 



da 



*A, 7 (0i,0 2 )d7 



(154) 



Lemma 5: Using the telescoping representation for x\(9), a solution for £a(0) is given as follows: 



(155) 



x A (6>) = / $\ !ti (9 1 a)x IJ ,(a)da+ / / $A )7 (0,a)du; 7 (a)da 
Je ./e 

^*\A d > «) = ^ [*a,„(0, a)] , $ A ,A = 5(0 - a) . 



Proof: We show that ( 11551 ) satisfies the differential equation in (f50b - Taking derivative of ( 11551 ) with 
respect to A, we have 



dx x (9) = J ^$a, m (0, a)x^{a)dad\ + J $a,a(#, a)du; A (a) + j_ j Jx $A ' 7 ^' a)dw 7 (a)da 



(156) 



/ F e [$ AiM (0,a)]^(a)(iadA + du;A(0)+ / / F e [$ A , 7 (0, a)] dw 7 (a)da (157) 



/e 



5> A)A( (0, a)x^(a)da + / / 3>> )7 (0, a)d6 7 (a)da 
Je Jn . 



F e [x x (e)} + dw x (e). 



+ dw x (6) 



(158) 
(159) 



To get (11581 ). we use Theorem [4] to take the integral transform Fq outside the integral. Since the x\(0) 
in dl55l ) satisfies (|50l ). it must be a solution. ■ 
Step 3. [Equation for x x \\(0)] Using (11381 ) and (11551 ). we can write g\ :fM {9,a) in (11471 ) as 



E x x {9) 



^7 



(; " <a) ^[x A (0)x 7 | 7 (a)]d 7 



d_ 

d_ 
d/j 
d_ 

dfi {Jo D 7 (a) ./„ 



G 7 (a). 



"" y -x 7 | 7 (a)d7 + n M (a) 



L> 7 (a) 



£>7(«) 

n G 7 (a) 



$ A , 7 (0, 0)£[s 7 (/3)5 7 i 7 (a)]d/3d7 



G M (a) / 



(160) 
(161) 
(162) 
(163) 
(164) 
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To get (11641 ). we use the fact that x^^a) _L x^^({3), so that 
Substituting (11641 ) in the expression for x X i x (6) in (11471) (Step 2), we get 



£aia(0) 



JB 



DJoc) 



e 



$ x ^9,f3)S^/3,a)d/3 



defj,(a)da . 



(165) 



Step 4. [Differential Equation for x A i A (6>)] Differentiating (11651 ) with respect to A, we get 



dx x \\(0) 



e D\(ot) 
Gx(a) 



S^O, a)de\(a)da + 



G^a) 



o Je D^a) Ue 
S\(a, 9)de x (a)da + Fg[x x \ x (9)]dX 



F e l$ x ^(9,l3)]S^f3,a)df3 



e -Da (« 
F e [x xlx (9)]dX + Kg[de x (d)}, 



de^(a)da 
(166) 
(167) 



where Kg is the integral transform defined as in (1581) . 

Step 5. [Differential Equation for S x (a,9)] The error covariance S x (a,9) can be written as 

S x (a, 9) = E[x xlx (a)x x{x (9)] = E[x x (a)x x (9)} - E[x x{x (a)x xlx (9)} . (168) 

Using expressions for x x (a) in ( 11551 ), we can show that for P x (a, f3) = E[x x (a)x x (9)], 

dP x (a,9) 



dX 



F a [P x (a, 9)} + F e P x (a, 9) + C x {a, 9) . 



(169) 



Using the expression for x x \ x {a) in (11461 ). it can be shown that 

dE[x x \ x (a)x x \ x {9)\ 



dX 



F a [ J B[x A | A (a)x A | A (0)]]+F e [ J B[x A | A (a)x A | A W]]+ / S x (a,P)S x (9, P)d? . 

(170) 



Differentiating (11681 ) and using (11691 ) and dl70| ), we get the desired equation 

GW) 



§^S x (a,9) = F a [S x (a,9)]+F [S x (a,9)}+C x (9,a) - J 



e DX/3) 



S x (a,f3)S x (9,P)dp. (171) 



Appendix D 
Proof of Theorem [9j Recursive Smoother 

We now derive smoothing equations. Using similar steps as in Lemma [3l we can show that 

-i 



x\\t(Q) 



g x j9,a)deja)da , 



(172) 



o Je 



DRAFT 



37 

gx, t ,(e,a) = -^E[x x (9)e t t(a)}. (173) 

Define the error covariance S\ :fl (9,a) as 

S^{6,a) = E[x x \ x {6)x^{a)]. (174) 

We have the following result for the smoother: 

Lemma 6: The smoothed estimator x\\t(9) is given by 

X\\t(Q) = + / / 9xA > a ) de n( a ) da > (175) 

Jx Je 



where for /x > A, 



gxAV^) = ^^S x ,,(d,a). (176) 



Proof: Equation ( 11751 ) immediately follows from d 1721 ) and Lemma [3j Equation d 1761) follows by 
using (11381 ) to compute gx,n(9,a) in (1 173b - ■ 
We now want to characterize the error covariance S x ^{9, a). Subtracting the telescoping representation 
in (150T ) and the filtering equation in d56l ), we get the following equation for the filtering error covariance: 

dx^a) = Felx^afldfi + B^{a)dw^{a) - K a [dn^(a)} , (177) 

where F a is the integral transform 



F a [xn\n{a)} = F a [x„\Ja)] - K a 



Just like we did in Lemma \5\ we can write a solution to (11771 ) as 
/<-) JeJx 



(178) 



z^{<*)= [ $ litX {a,6)x xlx (6)d6+ [ /"%, 7 (a, 0)[£ 7 (0)^ 7 (0) - K g [dn y {6)]]d6 (179) 
Je Je Jx 

^^x{a,e) = Fj>^x{a,0) and 0) = S(a - 9) . (180) 

Differentiating d 1791) with respect to A, we can show that 

/e #A j e 
Substituting (11791 ) in (11771) . we have the following relationship 



—$^ x (a,f3)x xlx ([3)df3 = - I <S>^ x (a,(3)F p x x \ x ((3)dp. (181) 



SxA°,<x)= / $n,\(a,P)Sx(P,e)dp. (182) 



e 
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Substituting (11821 ) in (11751 ) and (11761 ). differentiating (11751 ) and using (11811 ) and d59]), we get the 
following equation: 

dx MT (9) = Fg[x xlT (9)}dX+ f [ / ^ f ,,x(a,(3)C x ((3,e)d^de ll (a)da. (183) 

Jx Je L> u {a) j e 

Assuming S\(9,8) > 0, we get smoother equations using the following calculations: 

d£ X \ T (P)6(0-l3) = Ff i [x x \ T (p)]W-P)d\+ [ [ ^^r^„,x(a,P)C x (^e)de„(a)da 

Jx Je D u {a) 

(184) 

+ / / ^H^A^^xW,d)de^a)da (185) 

Sx(9,9) Sx(9,9) ,m U , 

lM9j) dx ^ {6) = C^9j) Fe[x ^ mdX 

+ J J^^^n^ lt ,x(a,P)S x iP,9)dp)de^a)da (186) 

dxx\ T {9) = F [xx\T(O)]d\ + ^^ I' [ ^\sxA0^)de,(a)da (187) 

S\{V,V) Jx Je D u (a) 

dx x{T (9) = F e [x x{T (9)]d\ + ^S^tW " ^A|aW] • (188) 

Equation ( fl83l is equivalent to (f1~84l. We multiply (PT841 by S^, 0) /C A (/3, 0) to get (PT85T ). We integrate 
(fl85T) for all /? to get (11861 To go from (fl86l) to (fT87l) . we use (fT82l) . Equation (TT87T) follows from (fr75l 
To derive a differential equation for S^^a, 0), we first note that 

S x]T (a,9) = E[x MT (a)x MT (9)} (189) 
= E[xx(a)xx(9)} - E[x x]T (a)x x]T (9)] . (190) 

Using ( 11751 ) to compute E[x x \t{o£)x x \t(0)], we can find an expression for S\| T (a,0) as 

Z" 1 f G 2 Ja) 

Sx\T(a,9) = S x {a,9) - / r- 5 M)A (ai, a)S^x(ai,9)dfida 1 . (191) 

Taking derivative of (1191b . we can derive (l66l) . 
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